%% PRACTICAL: OAT BEAMFORMING
%
% In this practical we will estimate neural activity in the brain's source
% space using a beamformer algorithm, task dependant differences in this source
% activity is then quantified using a GLM. This will go through the following
% steps:
% 
%   1) Set-up an OAT analysis: source_recon and first_level
%   2) Run source space GLM fitting and contrasts
%   3) Check coregistration
%   4) Save and view t-stat volumes
%   5) GLM analysis in a ROI
%   6) Time-frequency analysis in an ROI
%   7) Whole brain time-frequency analysis
%
% This practical will work with a single subject's data from an emotional
% faces experiment (Elekta Neuromag data). You can get the data from:
% https://sites.google.com/site/ohbaosl/practicals/practical-data/emotional-face-processing-elekta-neuromag-data
% 
% Work your way through the script cell by cell using the supplied dataset.
% As well as following the instructions below, make sure that you read all
% of the comments (indicated by %), as these explain what each step is
% doing. Note that you can run a cell (marked by %%) using the ?Cell? drop
% down menu on the Matlab GUI.    

%% SETUP THE MATLAB PATHS
%
% Sets the Matlab paths to include OSL. Change these paths so that they 
% correspond to the setup on your computer. You will also need to ensure 
% that fieldtrip and spm are not in your matlab path (as they are included
% within OSL).%% SETUP THE MATLAB PATHS
% make sure that fieldtrip and spm are not in your matlab path

global OSLDIR;
    
% This cell sets the Matlab paths to include OSL. Change the osldir path so
% that it corresponds to the setup on your computer before running the cell. 
osldir = '/Users/andrew/Software/Matlab/osl/osl_full/osl1.5.0_beta';

addpath(osldir);
osl_startup(osldir);

%% INITIALISE GLOBAL SETTINGS FOR THIS ANALYSIS
%
% This cell sets the directory that OAT will work in. Change the workingdir
% variable to correspond to the correct directory on your computer before
% running the cell.

% directory where the data is:
workingdir = '/Users/andrew/Projects/OSL_test/preproc_manual/faces_subject1_data';

cmd = ['mkdir ' workingdir]; if ~exist(workingdir, 'dir'), unix(cmd); end % make dir to put the results in

clear spm_files_continuous spm_files_epoched structural_files;

%% SET UP THE LIST OF SUBJECTS FOR THE ANALYSIS
%
% Specify a list of the fif files, structual files (not applicable for this
% practical) and SPM files (which will be created). It is important to make
% sure that the order of these lists is consistent across sessions. 
% Note that here we only have 1 subject, but more generally there would be
% more than one, e.g.:
%
% fif_files{1}=[testdir '/fifs/sub1_face_sss.fif']; 
% fif_files{2}=[testdir '/fifs/sub2_face_sss.fif']; 
% etc...
%
% spm_files{1} = [workingdir '/sub1_face_sss.mat'];
% spm_files{2} = [workingdir '/sub2_face_sss.mat'];
% etc...
% set up a list of SPM MEEG object file names (we only have one here)
% set up a list of SPM MEEG object file names (we only have one here)
% set up a list of SPM MEEG object file names (we only have one here)
%
% As we will be working in source space for this practical we will also define
% a list of structural MRI files, one for each participant. These will be used
% during coregistration and source localisation to define the source space for
% each participant
spm_files_continuous{1}=[workingdir '/spm8_meg1.mat'];
spm_files_epoched{1}=[workingdir '/espm8_meg1.mat'];

% set up a list of mris
structural_files{1}=[workingdir '/structurals/struct1.nii'];

%% SET UP WHOLE BRAIN BEAMFORMER IN OAT
%
% In previous practicals, the oat.source_recon parameters defined some simple
% data preparation to be done prior to sensor space analysis. In this practical
% we will define options for a whole brain source reconstruction  which
% requires some additional options.
%
% The conditions, freq range and time range are the same as we have used in
% previous analyses. The following parameters are used for the source analysis:
%   'method': This defines the source reconstruction method to be used. other options include 'beamform_bilateral' and 'mne_datacov'
%   'normalise_method': This defines how the data will be normalised
%   'gridstep': This is the distance (in mm) between points to be reconstructed, the spatial resolution of the analysis
%   'forward_meg': This is the forward model to be used
%   'mri': The path to the structural volume
%   'dirname': The output path


% SETUP THE OAT:
oat=[];
oat.source_recon.D_epoched=spm_files_epoched;
oat.source_recon.conditions={'Motorbike','Neutral face','Happy face','Fearful face'};
oat.source_recon.freq_range=[1 40]; % frequency range in Hz
oat.source_recon.time_range=[-0.2 0.4]; % time range in secs

oat.source_recon.method='beamform'; 
oat.source_recon.normalise_method='mean_eig';
oat.source_recon.gridstep=8; % in mm, using a lower resolution here than you would normally, for computational speed
oat.source_recon.forward_meg='Single Shell';

oat.source_recon.mri=structural_files;
oat.source_recon.dirname=[spm_files_continuous{1} '_erf_wideband_' oat.source_recon.method]; % directory the oat and results will be stored in

oat.source_recon.work_in_pca_subspace=0;


%% SETUP SOURCE SPACE OAT FIRST LEVEL GLM
%
% This cell defines the GLM parameters for the first level analysis, a
% trial-by-trial GLM analysis in the source space defined by oat.source_recon
% Critically this includes the design matrix (in Xsummary) and contrast
% matrix.
%
% Xsummary is a parsimonious description of the design matrix.
% It contains values Xsummary{reg,cond}, where reg is a regressor no. and cond
% is a condition no. This will be used (by expanding the conditions over
% trials) to create the (num_regressors x num_trials) design matrix:
%
% Each contrast is a vector containing a weight per condition defining how
% the condition parameter estimates are to be compared. Each vector will
% produce a different t-map across the sensors. Contrasts 1 and 2 describe
% positive correlations between the each sensors activity and the presence
% of a motorbike or face stimulus respectively. Contrast 3 tests whether
% each sensors activity is larger for faces than motorbikes.

Xsummary={};
Xsummary{1}=[1 0 0 0];
Xsummary{2}=[0 1 0 0];
Xsummary{3}=[0 0 1 0];
Xsummary{4}=[0 0 0 1];
oat.first_level.design_matrix_summary=Xsummary;

% contrasts to be calculated:
oat.first_level.contrast={};
oat.first_level.contrast{1}=[3 0 0 0]'; % motorbikes
oat.first_level.contrast{2}=[0 1 1 1]'; % faces
oat.first_level.contrast{3}=[-3 1 1 1]'; % faces-motorbikes
oat.first_level.contrast_name={};
oat.first_level.contrast_name{1}='motorbikes';
oat.first_level.contrast_name{2}='faces';
oat.first_level.contrast_name{3}='faces-motorbikes';
oat.first_level.report.first_level_cons_to_do=[2 1 3];
oat.first_level.time_range=[-0.1 0.3];
oat.first_level.post_tf_downsample_factor=1;
oat.first_level.name=['wholebrain_first_level'];

%% SANITY-CHECK OAT SETTINGS
%
% As well as using the settings we specified in the previous cell, calling osl_check_oat  has filled in a bunch of other settings as well.

oat = osl_check_oat(oat);

oat
oat.source_recon
oat.first_level


%% RUN THE OAT:
%
% The OAT structure you have created, oat, should be passed to the function
% osl_run_oat.m, to run the pipeline.  However, the OAT structure should
% contain a structure (oat.to_do), which is a list of binary values indicating
% which part of the pipeline is to be run. % E.g. oat.to_do=[1 1 0 0]; will run
% just the source recon and first level stages, whereas oat.to_do=[1 1 1 1];
% will run all four (source-recon, first level, subject-level and group level).
% Here we are not doing any group analysis, so we need:  oat.to_do=[1 1 0 0];

oat.to_do=[1 1 0 0];
oat = osl_run_oat(oat);

%% DO REGISTRATION AND RUN FORWARD MODEL BASED ON STRUCTURAL SCANS
%
% This step reruns the coregistration stage of the analyses.  This can just be
% done inside the oat.source_recon stage However, it is worth running
% separately to check that the results look reasonable!

S2=oat.source_recon;
S2.D = oat.source_recon.D_epoched{1}; 
S2.mri=oat.source_recon.mri{1}; 

D=osl_forward_model(S2);

spm_eeg_inv_checkdatareg(D);

%% OUTPUT SUBJECT'S NIFTII FILES
%
% Having run the GLM on our source space data, we would like to inspect the
% results for our single subject.  We can do this by saving the contrast of
% parameter estimates (COPEs) and t-statistics for each of our contrasts to
% NIFTI images.

S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.first_level_contrasts=[3]; % list of first level contrasts to output
S2.resamp_gridstep=oat.source_recon.gridstep;
[statsdir,times,count]=osl_save_nii_stats(S2);

%% VIEW NIFTII RESULTS IN FSLVIEW
%
% We can now view the nifti images containing our GLM results in FSL, here
% we are running fslview from the matlab command line, but you do not need
% to - you can run it from the UNIX command line instead.

mni_brain=[OSLDIR '/std_masks/MNI152_T1_' num2str(S2.resamp_gridstep) 'mm_brain']; 

% Inspect the results of a contrast in fslview
contrast=3;
runcmd(['fslview ' mni_brain ' ' [statsdir '/cope' num2str(contrast) '_' num2str(S2.resamp_gridstep) 'mm'] ' ' [statsdir '/tstat' num2str(contrast) '_' num2str(S2.resamp_gridstep) 'mm']  ' ' [statsdir '/tstat' num2str(contrast) '_mip_' num2str(S2.resamp_gridstep) 'mm'] ' &']);


%% INVESTIGATE A SINGLE LOCATION OF INTEREST 
% 
% In this section we will interrogate the wholebrain OAT (run above) using a
% specified MNI coordinate. This is sometimes referred to as a 'virtual
% electrode' or 'virtual sensor'
%
% The coordinate is defined in MNI152 space before the results for that
% location extracted from the OAT analysis we ran above

mni_coord=[27,-64,-18]; % set this to the MNI coord of interest

S2=[];
S2.vox_coord=mni_coord;
S2.stats=oat.first_level.results_fnames{1};
S2.oat=oat;
S2.first_level_cons_to_do=oat.first_level.report.first_level_cons_to_do; % plots all of these contrasts

[vox_ind_used] = oat_plot_vox_stats(S2);

%% INVESTIATE A REGION OF INTEREST
%
% In this section we will interrogate the wholebrain OAT (run above) using 
% an ROI mask
%
% Unlike the virtual electrode, the results from many voxels (all defined in
% the binary mask in S2.mask_fname) are extracted and the average results
% across these points presented.

% Spatially average the results over an ROI
S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.mask_fname=[OSLDIR '/std_masks/Right_Temporal_Occipital_Fusiform_Cortex'];
[stats,times,mni_coords_used]=oat_output_roi_stats(S2);

% plot
S2=[];
S2.stats=stats;
S2.oat=oat;
S2.first_level_cons_to_do=oat.first_level.report.first_level_cons_to_do; % plots all of these contrasts
[vox_ind_used] = oat_plot_vox_stats(S2);

%% ROI TIME-fREQUENCY ANALYSIS
%
% This section will re-rerun the first level analysis using a ROI in 
% the temporal occiptial fusiform cortex, and do a time-frequency
% trial-wise GLM first-level analysis
%
% We are going to use the wholebrain OAT (which was run above), to make use of the settings 
% and source_recon results already in there. The new parameters are
%   'tf_freq_range': the lower and upper bounds on the frequency range of interest
%   'tf_num_freqs': the number of frequency bands to estimate within the bounds in tf_freq_range
%   'tf_method': The spectral power estimation method
%   'tf_hilbert_freq_res': The resolution to use in the hilbert spectral estimation
%   'post_tf_downsample_factor': How much to downsample the tf results
%
% Give the first level analysis a new name to avoid copying over previous
% first-level analyses.

oat.first_level.name='roi_tf_first_level'; 
oat.first_level.mask_fname=[OSLDIR '/std_masks/Right_Temporal_Occipital_Fusiform_Cortex'];

oat.first_level.tf_freq_range=[4 48]; % frequency range in Hz
oat.first_level.time_range=[-0.1 0.3];
oat.first_level.tf_num_freqs=14;
oat.first_level.tf_method='hilbert';
oat.first_level.tf_hilbert_freq_res=8;
oat.first_level.post_tf_downsample_factor=2;

oat = osl_check_oat(oat);

oat.to_do=[0 1 0 0]; % We don't need to re-run the source_recon
oat = osl_run_oat(oat);

% VIEW RESULTS IN ANATOMICAL ROI
%
% The results within an anatomically defined right Fusiform region are then
% extracted and plotted

S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.mask_fname=[OSLDIR '/std_masks/Right_Temporal_Occipital_Fusiform_Cortex'];
[stats,times,mni_coords_used]=oat_output_roi_stats(S2);

% VIEW TF RESULTS FOR CONTRAST 3 (FACES>MOTORBIKES)
%
% Plot the time-frequency image for faces-motorbikes

S2=[];
S2.stats=stats;
S2.oat=oat;
S2.first_level_cons_to_do=3; % plots all of these contrasts
[vox_ind_used] = oat_plot_vox_stats(S2);

% VIEW RESULTS OVER TIME FOR A SINGLE FREQUENCY
%
% Plot the time course of a single freq bin

freqbin=nearest(stats.frequencies,8); % find bin for 8Hz

S2=[];
S2.stats=stats;
S2.oat=oat;
S2.freq_inds=freqbin;
S2.first_level_cons_to_do=3; % plots all of these contrasts
[vox_ind_used] = oat_plot_vox_stats(S2);

%% BROADBAND WHOLEBRAIN TF AMPLITUDE ANALYSIS
%
% This section will re-rerun the first level analysis over the whole brain, and
% do a time-frequency trial-wise GLM first-level analysis across a wide 4-40Hz
% band
%
% We are going to use the wholebrain OAT (which was run above for the ERF
% analysis), to make use of the settings already in there. This can be loaded
% with 'osl_load_oat_results'. The first-level time-frequency parameters are
% then modified

oat.source_recon.dirname=[spm_files_continuous{1} '_erf_wideband_beamform']; % directory the oat and results will be stored in
oat.first_level.name=['wholebrain_first_level'];
oat=osl_load_oat(oat);
try, oat.first_level=rmfield(oat.first_level,'freq_average'); catch, end;
res=osl_load_oat_results(oat,oat.source_recon.results_fnames{1});

% RUN THE OAT:
oat.to_do=[1 1 0 0];
oat.source_recon.gridstep=10; % in mm, using a lower resolution here than you would normally, for computational speed
oat.source_recon.freq_range=[4 40]; % broadband 

oat.first_level.tf_num_freqs=1;
oat.first_level.tf_method='hilbert';
oat.first_level.tf_freq_range=oat.source_recon.freq_range;
oat.first_level.post_tf_downsample_factor=2;
oat.first_level.name=['wholebrain_tf_first_level'];

oat = osl_check_oat(oat);
oat = osl_run_oat(oat);

%% VIEW RESULTS OF WHOLE BRAIN TF ANALYSIS
%
% Additional NIFTI volumes with different parameters (ie. gridspep or
% timepoints) can be saved out after the OAT analysis with the following code.
% This can then be checked in fslview

S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.first_level_contrasts=[3]; % list of first level contrasts to output
S2.freq_bin=1;
S2.stats_dir=[oat.source_recon.dirname '/' oat.first_level.name '_f' num2str(S2.freq_bin) '_stats_dir']; % directory to put niis in
S2.resamp_gridstep=oat.source_recon.gridstep;
[statsdir,times]=osl_save_nii_stats(S2);    

% view results using fslview
mni_brain=[OSLDIR '/std_masks/MNI152_T1_' num2str(S2.resamp_gridstep) 'mm_brain']; 
contrast_num=3;
runcmd(['fslview ' mni_brain ' ' statsdir '/tstat' num2str(contrast_num) '_' num2str(S2.resamp_gridstep) 'mm &']);
